Description: Load all required packages for data retrieval, manipulation, and plotting.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lubridate)
library(httr)
library(jsonlite)##
## Attaching package: 'jsonlite'
##
## The following object is masked from 'package:purrr':
##
## flattenlibrary(ggplot2)
library(purrr)
library(zoo) ##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(plotly) ##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:httr':
##
## config
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layoutDescription: Define the Trident bounding box and time window, then query the USGS FDSN API for all earthquakes in that region and time period.
# 1. Define Trident region and time window
trident_bbox <- list(
minlatitude = 58.0,
maxlatitude = 58.5,
minlongitude = -155.4,
maxlongitude = -154.8
)
starttime <- "2010-01-01"
endtime <- "2025-11-22" # update as needed
# 2. Build and send request to USGS
base_url <- "https://earthquake.usgs.gov/fdsnws/event/1/query"
resp <- GET(
url = base_url,
query = c(
format = "geojson",
starttime = starttime,
endtime = endtime,
trident_bbox,
eventtype = "earthquake",
minmagnitude = 0
)
)
stop_for_status(resp)
raw_json <- content(resp, "text", encoding = "UTF-8")
eq <- fromJSON(raw_json, flatten = TRUE)Description: Convert the raw GeoJSON features into a tidy event-level table with time, magnitude, location, and depth (in km and miles).
quakes <- eq$features %>%
as_tibble() %>%
transmute(
time = as.POSIXct(properties.time / 1000,
origin = "1970-01-01", tz = "UTC"),
mag = properties.mag,
place = properties.place,
lon = purrr::map_dbl(geometry.coordinates, 1),
lat = purrr::map_dbl(geometry.coordinates, 2),
depth_km = purrr::map_dbl(geometry.coordinates, 3), # depth from USGS
depth_miles = depth_km * 0.621371 # depth for public-facing labels
) %>%
arrange(time)
glimpse(quakes)## Rows: 10,408
## Columns: 7
## $ time <dttm> 2010-01-01 11:28:23, 2010-01-02 07:34:33, 2010-01-02 18:4…
## $ mag <dbl> 1.0, 1.0, 1.2, 1.6, 1.2, 1.7, 1.4, 1.6, 2.2, 1.2, 1.0, 1.6…
## $ place <chr> "99 km NNW of Karluk, Alaska", "86 km NW of Karluk, Alaska…
## $ lon <dbl> -154.9000, -155.3670, -154.9041, -155.3792, -154.9891, -15…
## $ lat <dbl> 58.4336, 58.1697, 58.4343, 58.1675, 58.2845, 58.1774, 58.1…
## $ depth_km <dbl> 1.0, 5.3, 1.0, 1.0, 2.7, 4.6, 2.8, 6.5, 1.0, 3.0, 3.2, 1.0…
## $ depth_miles <dbl> 0.6213710, 3.2932663, 0.6213710, 0.6213710, 1.6777017, 2.8…Description: Summarize event-level data into daily metrics that will feed into the unrest index (daily event count, magnitude distribution, and shallow-earthquake percentage). Negative depths are set to 0 km, and “shallow” is defined as < 3 km.
quakes_daily <- quakes %>%
mutate(
depth_km = ifelse(depth_km < 0, 0, depth_km), # clamp negative depths to 0 km
date = as.Date(time)
) %>%
group_by(date) %>%
summarise(
n_events = n(), # daily number of earthquakes
median_mag = median(mag, na.rm = TRUE), # typical quake size that day
p90_mag = quantile(mag, 0.9, na.rm = TRUE), # larger quake size indicator
shallow_pct = mean(depth_km <= 3, na.rm = TRUE), # % of shallow quakes (< 3 km)
.groups = "drop" # return ungrouped dataframe
)
head(quakes_daily)## # A tibble: 6 × 5
## date n_events median_mag p90_mag shallow_pct
## <date> <int> <dbl> <dbl> <dbl>
## 1 2010-01-01 1 1 1 1
## 2 2010-01-02 3 1.2 1.52 0.667
## 3 2010-01-09 1 1.2 1.2 1
## 4 2010-01-12 1 1.7 1.7 0
## 5 2010-01-14 1 1.4 1.4 1
## 6 2010-01-15 2 1.9 2.14 0.5Description: Configure how the baseline is defined. • baseline_mode = “fixed” uses a pre-unrest historical window. • baseline_mode = “rolling” uses a rolling window (previous N years) for each date. • baseline_mode = “anchored” uses a 3/5/7-year or historical window before a chosen reference date t₀.
# -------------------------------------------------------------------
# Baseline configuration
# -------------------------------------------------------------------
# Choose baseline mode: "fixed", "rolling", or "anchored"
baseline_mode <- "anchored"
# For rolling baseline (not the main focus right now)
baseline_window_years <- 3 # previous N years per date if baseline_mode == "rolling"
# Fixed (historical) baseline: user-defined quiet period (pre-unrest)
baseline_fixed_start <- as.Date("2015-01-01")
baseline_fixed_end <- as.Date("2022-08-24") # start of documented unrest at Trident
# Anchored baseline: history before a chosen reference date t0
anchored_ref_date <- as.Date("2022-08-24") # t0; can be any date of interest
anchored_span_years <- 5 # 3, 5, or 7 years
anchored_historical <- FALSE # TRUE = full history up to t0 (ignore span)Description: For rolling mode, compute per-date baseline statistics using the previous N years of daily seismic metrics (excluding the current day).
compute_rolling_baseline <- function(quakes_daily,
window_years = 3,
min_days = 30) {
dates <- quakes_daily$date
purrr::map_dfr(dates, function(d) {
start <- d %m-% years(window_years) # window_years before date d
ref <- quakes_daily %>%
filter(date >= start, date < d) # strictly before current day
if (nrow(ref) < min_days) {
tibble(
date = d,
n_mu = NA_real_,
n_sd = NA_real_,
p90_mu = NA_real_,
p90_sd = NA_real_,
shp_mu = NA_real_,
shp_sd = NA_real_
)
} else {
tibble(
date = d,
n_mu = mean(ref$n_events, na.rm = TRUE),
n_sd = sd(ref$n_events, na.rm = TRUE),
p90_mu = mean(ref$p90_mag, na.rm = TRUE),
p90_sd = sd(ref$p90_mag, na.rm = TRUE),
shp_mu = mean(ref$shallow_pct, na.rm = TRUE),
shp_sd = sd(ref$shallow_pct, na.rm = TRUE)
)
}
})
}Description: • Standardize each daily metric relative to its baseline (fixed, rolling, or anchored) using z-scores. • Convert z-scores to intensities (0–1), average them, and rescale to 0–100. • Assign a traffic-light status (Green/Yellow/Orange/Red) based on the unrest score.
intensity_from_z <- function(z) {
x <- z / 3 # 3 sd ~= "max" intensity
x <- pmin(pmax(x, 0), 1)
return(x)
}
eps <- 1e-6 # to avoid divide-by-zero
if (baseline_mode == "fixed") {
# ----------------------------------------------
# Fixed baseline: single historical quiet period
# ----------------------------------------------
baseline_stats <- quakes_daily %>%
filter(date >= baseline_fixed_start,
date < baseline_fixed_end) %>%
summarise(
n_mu = mean(n_events, na.rm = TRUE),
n_sd = sd(n_events, na.rm = TRUE),
p90_mu = mean(p90_mag, na.rm = TRUE),
p90_sd = sd(p90_mag, na.rm = TRUE),
shp_mu = mean(shallow_pct, na.rm = TRUE),
shp_sd = sd(shallow_pct, na.rm = TRUE)
)
quakes_scored <- quakes_daily %>%
mutate(
z_n = (n_events - baseline_stats$n_mu) / (baseline_stats$n_sd + eps),
z_p90 = (p90_mag - baseline_stats$p90_mu) / (baseline_stats$p90_sd + eps),
z_shp = (shallow_pct - baseline_stats$shp_mu) / (baseline_stats$shp_sd + eps),
I_n = intensity_from_z(z_n),
I_p90 = intensity_from_z(z_p90),
I_shp = intensity_from_z(z_shp),
unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
status = case_when(
unrest_score < 20 ~ "Green",
unrest_score < 40 ~ "Yellow",
unrest_score < 70 ~ "Orange",
TRUE ~ "Red"
)
)
} else if (baseline_mode == "rolling") {
# ----------------------------------------------
# Rolling baseline: previous N years per date
# ----------------------------------------------
baseline_rolling <- compute_rolling_baseline(
quakes_daily,
window_years = baseline_window_years
)
quakes_scored <- quakes_daily %>%
left_join(baseline_rolling, by = "date") %>%
mutate(
z_n = (n_events - n_mu) / (n_sd + eps),
z_p90 = (p90_mag - p90_mu) / (p90_sd + eps),
z_shp = (shallow_pct - shp_mu) / (shp_sd + eps),
I_n = intensity_from_z(z_n),
I_p90 = intensity_from_z(z_p90),
I_shp = intensity_from_z(z_shp),
unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
status = case_when(
unrest_score < 20 ~ "Green",
unrest_score < 40 ~ "Yellow",
unrest_score < 70 ~ "Orange",
TRUE ~ "Red"
)
)
} else if (baseline_mode == "anchored") {
# ----------------------------------------------
# Anchored baseline: history before a reference date t0
# ----------------------------------------------
# Determine baseline window relative to anchored_ref_date
if (anchored_historical) {
baseline_start <- min(quakes_daily$date, na.rm = TRUE)
} else {
baseline_start <- anchored_ref_date %m-% years(anchored_span_years)
}
baseline_end <- anchored_ref_date
baseline_stats <- quakes_daily %>%
filter(date >= baseline_start,
date < baseline_end) %>%
summarise(
n_mu = mean(n_events, na.rm = TRUE),
n_sd = sd(n_events, na.rm = TRUE),
p90_mu = mean(p90_mag, na.rm = TRUE),
p90_sd = sd(p90_mag, na.rm = TRUE),
shp_mu = mean(shallow_pct, na.rm = TRUE),
shp_sd = sd(shallow_pct, na.rm = TRUE)
)
quakes_scored <- quakes_daily %>%
mutate(
z_n = (n_events - baseline_stats$n_mu) / (baseline_stats$n_sd + eps),
z_p90 = (p90_mag - baseline_stats$p90_mu) / (baseline_stats$p90_sd + eps),
z_shp = (shallow_pct - baseline_stats$shp_mu) / (baseline_stats$shp_sd + eps),
I_n = intensity_from_z(z_n),
I_p90 = intensity_from_z(z_p90),
I_shp = intensity_from_z(z_shp),
unrest_score = 100 * (I_n + I_p90 + I_shp) / 3,
status = case_when(
unrest_score < 20 ~ "Green",
unrest_score < 40 ~ "Yellow",
unrest_score < 70 ~ "Orange",
TRUE ~ "Red"
)
)
}Description: To reduce noise from day-to-day fluctuations in seismicity, a 7-day rolling mean is applied to the unrest index. This smoothing highlights broader trends in volcanic activity, making the seismic state easier to interpret in dashboards and public-facing visualizations. The raw daily score is preserved for internal analysis, while the averaged score (unrest_score_7d) is used for plotting.
quakes_scored <- quakes_scored %>%
arrange(date) %>%
mutate(
unrest_score_7d = as.numeric(zoo::rollapply(
unrest_score,
width = 7,
FUN = mean,
align = "right",
fill = NA,
na.rm = TRUE
))
)
summary(quakes_scored$unrest_score_7d)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.9969 10.3291 13.8954 16.9065 20.1364 50.8806 6Description: Compute the mean of the 7-day smoothed unrest score over the chosen baseline window (3/5/7 years or historical before the reference date). This will be shown as a dashed horizontal line on the plot.
if (baseline_mode == "fixed") {
baseline_start_plot <- baseline_fixed_start
baseline_end_plot <- baseline_fixed_end
} else if (baseline_mode == "anchored") {
if (anchored_historical) {
baseline_start_plot <- min(quakes_scored$date, na.rm = TRUE)
} else {
baseline_start_plot <- anchored_ref_date %m-% years(anchored_span_years)
}
baseline_end_plot <- anchored_ref_date
} else {
baseline_start_plot <- NA
baseline_end_plot <- NA
}
if (!is.na(baseline_start_plot)) {
baseline_mean_score <- quakes_scored %>%
filter(date >= baseline_start_plot,
date < baseline_end_plot) %>%
summarise(mu = mean(unrest_score_7d, na.rm = TRUE)) %>%
pull(mu)
} else {
baseline_mean_score <- NA_real_
}
baseline_start_plot## [1] "2017-08-24"baseline_end_plot## [1] "2022-08-24"baseline_mean_score## [1] 12.90601Description: Visualize the 7-day smoothed unrest score over time with color-coded status, a vertical line marking the anchored reference date, and a dashed horizontal line showing the baseline mean over the chosen window.
subtitle_text <- dplyr::case_when(
baseline_mode == "fixed" ~ paste0(
"Fixed baseline: ", baseline_fixed_start, " to ", baseline_fixed_end, " (pre-unrest)"
),
baseline_mode == "rolling" ~ paste0(
"Rolling baseline: previous ", baseline_window_years, " years per day"
),
baseline_mode == "anchored" & anchored_historical ~ paste0(
"Anchored historical baseline up to ", anchored_ref_date
),
baseline_mode == "anchored" ~ paste0(
"Anchored baseline: ", anchored_span_years, "-year window before ", anchored_ref_date
),
TRUE ~ ""
)
# Vertical line date: baseline_fixed_end for fixed, anchored_ref_date for anchored
vline_date <- if (baseline_mode == "fixed") {
baseline_fixed_end
} else if (baseline_mode == "anchored") {
anchored_ref_date
} else {
NA
}
p <- ggplot(
quakes_scored,
aes(
x = date,
y = unrest_score_7d,
color = status,
text = paste0(
"Date: ", date,
"<br>Score (7d): ", round(unrest_score_7d, 1),
"<br>Daily events: ", n_events,
"<br>p90 mag: ", round(p90_mag, 2),
"<br>Shallow %: ", round(shallow_pct * 100, 1), "%"
)
)
) +
geom_line(na.rm = TRUE) +
geom_vline(
xintercept = as.numeric(vline_date),
linetype = "dashed"
) +
{
if (!is.na(baseline_mean_score)) {
geom_hline(yintercept = baseline_mean_score,
linetype = "dashed")
} else {
NULL
}
} +
# Baseline label annotation (rounded to 2 decimal places, placed above the line)
{
if (!is.na(baseline_mean_score)) {
annotate(
"text",
x = max(quakes_scored$date, na.rm = TRUE),
y = baseline_mean_score + 2, # 2 units above the line on the 0–100 scale
label = paste0("Baseline \u2248 ", round(baseline_mean_score, 2)),
hjust = -0.1,
vjust = 0.5,
size = 3.5,
color = "black"
)
} else {
NULL
}
} +
labs(
title = "Trident Volcano – 7-Day Seismic Unrest Score",
subtitle = subtitle_text,
x = "Date",
y = "Unrest Score (0–100, 7-day mean)",
color = "Status"
) +
scale_color_manual(
values = c(
"Green" = "#00CC00",
"Yellow" = "#FFD700",
"Orange" = "#FF9900",
"Red" = "#CC0000"
),
breaks = c("Red", "Orange", "Yellow", "Green")
) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "#F0F0F0", color = NA),
plot.background = element_rect(fill = "#F0F0F0", color = NA),
panel.grid.major = element_line(color = "grey80"),
panel.grid.minor = element_line(color = "grey90"),
clip = "off"
)Description: The final plot is rendered as an interactive HTML widget using plotly::ggplotly(). Hovering over the line (and points) reveals detailed values for each date, including the 7-day unrest score, daily earthquake count, magnitude distribution, and shallow earthquake percentage.
# Add points to make individual days easy to see in the interactive view
p <- p + geom_point(size = 0.5, alpha = 0.9)
ggplotly(p, tooltip = "text")## Warning in scale_x_date(): A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.## Warning in plot_theme(list(theme = theme), default = default): The `clip` theme
## element is not defined in the element hierarchy.